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O ' Abstract 

Von Neuman 's work on universal machines and the hardware development have allowed the simulation of dynamical 
, systems through a large set of interacting agents. This is a bottom-up approach which tries to derive global properties of a 
complex system through local interaction rules and agent behaviour Traditionally, such systems are modeled and simulated 
through top-down methods based on differential equations. Agent-Based Modeling has the advantage of simplicity and low 
computational cost. However, unlike differential equations, there is no standard way to express agent behaviour Besides, 
it is not clear how to analytically predict the results obtained by the simulation. Such observations got the attention of the 
scientific community and some techniques have been proposed in order to cover these gaps in the agent-based modeling field. 
In this paper we survey some of these methods. For expressing agent behaviour formal methods, like Stochastic Process 
Algebras have been used. Such approach is useful if the global properties of interest can be expressed as a function of 
. stochastic time series. However, if space variables must be considered, that means, if the space distribution of agents is 
C3 ' important we shall change the focus. In this case, multiscale techniques, based on Chapman-Enskog expansion was used to 
establish the connection between the microscopic dynamics ( agent behaviour) and the macroscopic observables. Besides, 
knowledge discovery in agent systems is a NP problem. This is the motivation for using data mining techniques, like Principal 
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. Component Analysis (PCA), to study agent systems like Cellular Automata. With the help of these tools (Stochastic Process 

' Algebras, Chapman-Enskog expansion and PCA) we will discuss a simple society model, a Lattice Gas Automaton for fluid 

^ modeling, and knowledge discovery in CA databases. Besides, we show the capabilities of the NetLogo, a free software for 

' agent simulation of complex system and describe our experience with this package. 
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c5 1. Introduction 



With the development of the hardware the possibility of simulating a system by constructing a mathematical model and 
executing it on a computer has opened new frontiers in science and engineer |22 31 10 30). Traditionally, the mathematical 
model is based on differential equations connecting the macroscopic variables that define the system 1221 13 II . For example, 
the majority of the fluid models follow the Eulerian formulation of fluid mechanics; that is, the fluid is considered as a 
continuous system subjected to Newton's and conservation Laws as well as state equations connecting the thermodynamic 
variables of pressure P, density p and temperature T I22I . This is a top-down approach which attempts to capture the nature 
of the relationships between macroscopic variables without been specific about the essence of the microscopic scales. 

On the other hand, agent-based modeling tries to emulate the system behavior following another viewpoint |2 30j. In 
this case, the model consists of a set of agents that encapsulate the behaviors of the individuals that make up the system, 
and execution consists of emulating these behaviors f7"3"21 1. These are bottom-up models based on the description of the 
individuals (agents) and their local interactions as well as the belief that the macroscopic observables and their relationships 
can be derived from the microscopic (agents) interactions. For instance, that is the philosophy behind Lattice Gas Cellular 
Automata models for fluids f 11 1 as well as some techniques for simulating social and ecological processes 0). 



In this paper we focus on agent-based models for natural phenomena. We observe two approaches in this field: Cellular 
Automata and Agent-Based Cellular Automata approaches. Cellular Automata are discrete and finite dynamical systems 
that evolve following simple and local rules which can be deterministic of probabilistic ones. For example, in modeling 
pheromone trails 1321 . each cell might contain a pair of state values as well as the amount of pheromone at a certain position 
and a binary value determining whether or not an ant is present in that cell. If a cell contains an ant then it will move to 
the adjoining cell with the most pheromone, depositing pheromone in the cell it leaves. Otherwise, the pheromone in a cell 
without an ant will decrease (due to evaporation). 

Instead of expressing the rules of the above model in terms of update rules for cells, the rules could be equally well 
expressed in terms of how each ant behaves, that is, an algorithm is used to describe the behavior of the ant and if it moves 
between cells, on each time step choosing the neighbouring cell with the most pheromone. In this viewpoint, the agent- 
base one, we can abstract the space distribution of agents and focus in their activities and interactions. Obviously, space 
distributions are easily recovered by imposing that ants move on a lattice. So agent based modeling incorporates the cellular 
automata philosophy also. 

Agent-Based Modeling, has the advantage of simplicity and low computational cost. However, unlike differential equa- 
tions, there is no standard way to express agent behavior Besides, it is not clear how to analytically predict the results 
obtained by the simulation. 

This paper is organized as follows. The next section presents the basic concepts of CAs and how computational intractable 
problems arise in this area. Then, Section |3l shows the application of PCA for cellular automata analysis. In Section|5]we 
review the WSCCS, a stochastic process algebra, and its application for expressing agent behavior and interaction. Section|S] 
presents the Chapman-Enskog expansion in the context of cellular automata for fluid modeling. In Section0we describe the 
NetLog capabilities and present our implementation of the HPP through NetLog tools. Finally, we discuss some perspectives 
in the field of agent-based modeling and simulation. 

2. Cellular Automata 

A cellular automaton (CA) is a quadruple A = {L; S; N; /) where L is a set of indices or sites, S is the finite set of site 
values or states, TV : L — > L*^ is a one-to-many mapping defining the neighborhood of every site i as a collection of k sites, 
and f : S'^ ^ S is the evolution function of A 1 38 5 1. The neighborhood of site i is defined as the set N{i) — {j; \ j — i\ < 
[{k — l)/2]} ([x] stands for the integer part of x). Note that a given site may or not be included in its own neighborhood. 
Since the set of states is finite, {fj} will denote the set of possible rules of the CA taken among thep = (#S')(#^)' rules. 

For a one-dimensional cellular automaton the lattice L is an array of sites, and the transition rule / updates a site value 
according to the values of a neighborhood of fc = 2r + 1 sites around it, that means: 

/ : S'2'-+i -> S, (1) 



a* e S, j ^i-r,...,i + r. (3) 

where t means the evolution time, also taking discrete values, and a' means the value of the site i at time t fW4] (see also 
(36 1 for on-line examples). Therefore, given a configuration of site values at time t, it will be updated through the application 
of the transition rule to generate the new configuration at time t+l, and so on. In the case of r = 1 in Expression (|2ji and 
S = {0, 1} we have a special class of cellular automata which was widely studied in the CA literature I29lll2l l9l l38l . Figure 
[T]shows the very known example of such a CA. The rule in this case is: 

= (a*_i + a*+i) mod2, (4) 

that means, the remainder of the division by two. The figure pictures the evolution of an initial configuration in which there 
is only one site with the value 1 . 

Once r = 1 in Expression (|2ji, it is easy to check that this rule is defined by the function: 



111 110 101 100011 010 001 000 
01011010 



Figure 1. Evolution of CA given by rule in expression!?] In this case, the initial configuration is a 
finite one-dimensional lattice which has only one site with the value 1 (pictured in black). 



x= 2^ + 1 * 2^ + * 2^ + 1 * 2"* + 1 * 2^ + * 2^ + 1 * 2^ + * 2° = 90 (6) 

By observing this example, we see that there are 2^ — 256 such rules and for each one it can be assigned a rule number 
following the indexation illustrated on Expression (|6}. In |37|, Wolfram proposes four basic classes of behavior for these 
rules (see also |5 1): 

Class 1: Evolution leads to homogeneous state in which all the sites have the same value (Figure|2]a); 
Class 2: Evolution leads to a set of stable and periodic structures that are separated and simple (Figure|2]b); 
Class 3: Evolution leads to a chaotic pattern (Figure|2lc); 
Class 3: Evolution leads to complex structures (Figure|2ld). 




(c) 




Figure 2. Some examples of Wolfram's classification for one-dimensional (r = 1) CAs. 

Other classifications based on Markovian processes and group properties can be also found in the literature I20II13I . 



Despite of its local simplicity, knowledge discovery in CA is a NP problem. In fact, let us take a one-dimensional CA with 
a finite lattice L of size d. One may consider the question of whether a particular sequence of d site values can occur after T 
time steps in the evolution of the cellular automaton, starting from any initial state. Then, one may ask whether there exists 
any algorithm that can determine the answer in a time given by some polynomial in d and T. The question can certainly be 
answered by testing all sequences of possible initial site values, that is {^SY- But this procedure requires a time that grows 
exponentially with d. 

Nevertheless, if an initial sequence could be guessed, then it could be tested in a time polynomial in d and T. As 
a consequence, the problem is in the class NP which motivates the application of data mining techniques for knowledge 
discovery in CA. The next sections review PCA basic theory and its application for the analysis of the (traditional) set of 
rules composed by 256 ID cellular automata obtained when r = 1, S* = {0, 1}. 

3. Principal Component Analysis 

Principal Component Analysis (PCA), also called Karhunen-Loeve, or KL method, can be seen as a method for data 
compression or dimensionality reduction f6^ (see f24l, section 5.11 also). Thus, let us suppose that the data to be compressed 
consist of N tuples or data vectors, from a n-dimensional space. Then, PCA searches for k n-dimensional orthonormal 
vectors that can best be used to represent the data, where k < n. Figure |3la-b pictures this idea using a bidimensional 
representation. If we suppose the data points are distributed over the ellipse, it follows that the coordinate system ((X, Y) 
shown in Figure|3]b) is more suitable for representing the data set in a sense that will be formally described next. 

Thus, let S = {ui, U2, Ujy} be the data set represented on Figure |3] By now, let us suppose that the centroid of the 
data set is the center of the coordinate system, that means: 

1 ^ 

Cm = ^5I"* = 0. (7) 




Figure 3. (a)Original dataset. (b) Extraction of the principal component. 

To address the issue of compression, we need a vector basis that satisfies a proper optimization criterion (rotated axes in 
Figure|3]b). Following |24|, consider the operations in Figure^] The vector Uj is first transformed to a vector by the 
matrix (transformation) A. Thus, we truncate by choosing the first m elements of Vj . The obtained vector Wj is just the 
transformation of Vj by /„, that is a matrix with Is along the first m diagonal elements and zeros elsewhere. Finally, Wj is 
transformed to zj by the matrix B. Let the square error defined as follows: 
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where Tr means the trace of the matrix between the square brackets and the notation ( *T) means the transpose of the complex 
conjugate of a matrix. Following Figure^] we observe that Zj = BImAuj. Thus we can rewrite (|8j as: 
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Figure 4. KL transform formulation. 
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which yields: 



where: 
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iI~BI^A)RiI-BIraA)*^ 
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(9) 
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i=0 



Following the literature, we call R the covariance matrix. We can now stating the optimization problem by saying that we 
want to find out the matrices A, B that minimizes J,„. The next theorem gives the solution for this problem. 
Theorem 1: The error J,„ in expression ilO\ is minimum when 



A = $ 



*T 



B = AB ^ BA^ L 



(12) 



where $ is the matrix obtained by the orthonormahzed eigenvectors of R arranged according to the decreasing order of its 
eigenvalues. 

Proof. To minimize J,,, we first observe that J„, must be zero if m = n. Thus, the only possibility would be 



Besides, by remembering that 



we can also write: 



I ^ BA^ A = B- 



Tr (CD) = Tr (DC) , 



(13) 



(14) 



T — —Tr 

n 

Again, this expression must be null \f m = n. Thus: 

1, 



{I~BI^A)*^ {I-BIraA)R 



(15) 



This error is minimum if: 



J„ = -Tr \(I - BA- A*^B*^ + A*^B*^BA) R] 
n 



B*^B = I, A*^A = I, 



(16) 



that is, if A and B are unitary matrix. The next condition comes from the differentiation of J„i respect to the elements of A. 
We should set the result to zero in order to obtain the necessary condition to minimize J„j. This yields: 



which renders: 



I^A*^ {I - A*^I,nA) R^O, 



J„, = -Tr [{I - A*^I,nA) R] 



(17) 
(18) 



By using the property (I14> . the last expression can be rewritten as 



Jra = -Tr [R - Ir^ARA*^] . 



Since R is fixed, J„, will be minimized if 



J™ = Tr [I^ARA*^] = J2 



(19) 



is maximized where af is the ith row of A. Once A is unitary, we must impose the constrain: 



T * 1 

a, a, = 1. 



(20) 



Thus, we shall maximize J,„ subjected to the last condition. The Lagrangian has the form: 



i=0 i=0 

where the are the Lagrangian multipliers. By differentiating this expression respect to a; we get: 

Rci* — Aid* , 

Thus, a* are orthonormalized eigenvectors of R. Substituting this result in expression ( I19I I produces: 

m — 1 



— ^ ^ Ai 



i=0 



(21) 



(22) 



which is maximized if {a*, i = 0,1, ...,m — 1} correspond to the largest eigenvalues of R. (□) 

A straightforward variation of the above statement is obtained if we have a random vector u with zero mean. In this case, 
the pipeline of Figure|3]yields a random vector z and the square error can be expressed as: 



J 17 



-Tr 



E^iu- BIraAu) (U - BlraAu)*^^ 



which can be written as: 



Jm = -Tr 
n 



{I-BI„,A)R{I-BI„,Ay 



(23) 



where R — E (uu*^) is the covariance matrix. Besides, if Cm in Expression is not zero, we must translate the coordinate 
system to Cm before computing the matrix R , that is: 



Uj = Uj-C„ 



(24) 



In this case, matrix R will be given by: 



N 



Also, sometimes may be useful to consider in Expression some other norm, not necessarily the 2-norm. In this case, there 
will be a real, symmetric and positive-defined matrix M, that defines the norm. Thus, the square error J,„ will be rewritten 
in more general form: 



1 ^ 2 1 ^ T 



(25) 



j=o 



3=0 



Obviously, if M = / we recover Expression (|8j. The link between this case and the above one is easily obtained by 
observing that there is non-singular and real matrix W , such that: 



W^MW = I. 



(26) 



The matrix W defines the transformation: 



Thus, by inserting these expressions in Equation (125 > we obtain: 



1 ^ ^ *T ^ 



Expression ( I28l l can be written as: 



1 ^ 



II > 



now using the 2-norm, like in Expression (|8}. Therefore: 



Jrn = -Tr 



N 



3=0 



Following the same development performed above, we will find that we must solve the equation: 

Ra* = Aia*, 



where: 



AT 



Thus, from transformations J27> it follows that: 

R = WRW'^. 

and, therefore, we must solve the following eigenvalue/eigenvector problem: 

{WRW^)a* = X,a*. 
The eigenvectors, in the original coordinate system, are finally given by: 

Wa* = a*. 

The next section shows the application of PCA method for knowledge discovery in CAs. 
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4. PCA and Cellular Automata 

In this section we review the work presented in [14]. In this reference, the authors analyzed one-dimensional CAs using 
PCA. The key idea is to consider binary patterns of a pre-defined size I as inputs of the CAs. It is considered the 256 one- 
dimensional CA rules obtained for r — 1 and S — {0, 1} in expression ! 1121 The output can be collected in a Table, like Table 
□ built for I = 5. 
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00000 


000 


Ill . 


. 000 


111 


00001 


000 


110 . 


. 001 


111 


11110 


000 


000 . 


. Ill 


111 


11111 


000 


000 . 


. Ill 


111 



Table 1. Table which rows are indexed by binary patterns and collumns by the CA rules Ro, Ri, 

R255- 



Each row j of Table ^ is obtained through the application of the rule Rj (see Expression (|6j for an example of rule 
indexation) Then, I/O patterns are converted to cardinal numbers denoted by fj {nii), which means the cardinal number 
corresponding to the application of the rule j to the pattern i {i — 0,1, ...,31 for Table^. Thus, in general, we get the matrix: 



F 



h 



In 



hp 
fnp 



(36) 



where fij = fj {rrii) . The matrix F is the data set to be analyzed. 

For mining knowledge in F through PCA we should firstly to perform the operation (translation) given by Thus, 
matrix F is converted to the following one: 



X = 



Xii 



with: 



(37) 



(38) 



E3 = -y2fjim.d. (39) 

i=l 

The matrix X is of size np. In 1141 columns xi, . . . ,Xp of X are called variables while rows ei , . . . , e„ are called covariables. 
However, we must observe that space dimension is the number of rules (p) and the number of data vectors is the number of 
patterns (n). Thus, following Section|3 we should apply the PCA over the data set given by matrix X'^ in order to find out 
the principal components of the covariables space. Besides, in 1141 the norm of the covariables space is defined by: 

M ^dtag(^-^,j^,...,j^y (40) 

with: 

S^-^i-t4-(t-.)] = '-ti-.-E,f- (41) 

\ i=l \i=l / I i=l 

Following Section |3] we must solve Equation ( I34> to find the eigenvalues and then apply Expression i35i to get the 
eigenvectors in the desired representation. The Table|2]shows the larger eigenvalues of this matrix for the listed pattern sizes. 

The main result is that the eigenvalues from the seventh rank are dramatically smaller in magnitude (104 times) than the 
first seven ones. Such observation led authors of 1 14| towards the following conjecture: 

Conjecture: The rank of R is 7 and does not depend on the size I of patterns being considered. When I is increased the 
eigenvalues tend to characteristic values obtained for I = 12. 

This is the main result presented in L14J . Next, we show our results by applying the same analysis but introducing 
randomness in the CA behavior. 
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Ai 


A2 


A3 


A4 


As 


Ae 


Ay 


4 


52.6802 


48.2214 


36.8869 


36.8263 


36.3134 


24.4539 


18.6179 


5 


58.2575 


50.9776 


37.2301 


37.0399 


30.7382 


21.7355 


18.0214 


6 


59.5952 


51.6519 


37.3406 


37.1109 


29.3769 


21.0940 


17.8305 


7 


59.9260 


51.8197 


37.3696 


37.1296 


29.0383 


20.9358 


17.7811 


9 


60.0290 


51.8721 


37.3788 


37.1355 


28.9325 


20.8865 


17.7656 


12 


60.0358 


51.8755 


37.3794 


37.1359 


28.9256 


20.8833 


17.7645 



Table 2. Eigenvalues of the correlation matrix. 
5. Stochastic Process Algebra and Agents 

If we can break down a system into component parts that act as finite state machines, then we can apply formal methods 
to explain how they are combined to form the observed (macroscopic) whole. That is the key idea of using process algebra 
for modeling societies 121ll30l l. Process algebra are widely used in the analysis of distributed computer systems (I]. They 
allow formal reasoning about how the various components of a system contribute to its overall behavior 1261 1231 . 

In 1 32 1, it is argued that a stochastic process algebras, the Weighted Synchronous Calculus of Communicating Systems 
(WSCCS), provides a useful formalism for understanding the dynamical behavior of their colony, since they combine com- 
puter simulation, Markov chain analysis and mean-field methods of analysis. Next, we review the basic elements of a process 
algebra and show its application for modeling societies. 

5.1. Stochastic Process Algebra 

One of the best known process algebra, and also a remarkable one in this area, is the Calculus of Communicating Systems 
(CCS) |1 1. It uses the notions of agents (or processes) and actions. Agents describe the entities which make up a system, 
such as processes in a distributed system, and actions that allow the agents communication (interaction). These notions are 
formally described which permits logical reasoning about the system iSj. Besides, in a new equivalence concept for 
agents, which are finite state automata, is provided. The CCS makes no attempt to actions synchronization and priority. The 
WSCCS adds such features to the CCS |33 34 35 1. 

Any process algebra consists of essentially four components [8 j: 

1 . A syntax for describing agents (automata) and the actions they perform. 

2. Algebraic rules. 

3. Derivation rules. 

4. A congruence for defining when two automata are considered equivalent 
in all algebraic contexts. 

5. An equational theory which defines how the equivalence of automata is demonstrated from the syntax of the agents 
which compose them. 

For instance, in the WSCCS it is used the following syntax: 
Agents are labeled by capital letters like A, B, C, ... 

The set of allowed actions Act form an abelian group {Act, *), where * is the group operation. The identity action, denoted 
by /, can be seen as a tick of a global clock. Each time a / occurs time has just moved forward one step. The inverse of an 
action a is denoted by a, which means, a*a= / . This operation will formally represent communication between agents in 
the WSCCS. 

For example, let us suppose that we have two agents A and B and that, in a single unit of time, there is a probability p that 
A becomes B and a probability 1 — p that it remains unchanged. Thus, we can define A by the following algebraic expression 
in WSCCS: 

A = j>:/.B + {\-v):^.A, (42) 

where the + indicates that the agent can make a choice. In Expression J42t . each possible choice will define a transition and 
the transitions will define the derivation rules. Formally, we write: 



A /[p] B. 



(43) 



In general, we have: 



A /[l-p] A. 

> 



(44) 



A a[p\ B (45) 

which means that agent A may change to B, with probability p, when action a occurs. 

Another important operation is the composition (x) of agents. Given the agents a.E and h.F, where a, h are possible 
actions, their composition is formally defined by: 

a.E X h.F = ah. {E x F) . (46) 
This expression do not incorporates the probability. The following expression adds this feature: 



(^a,:E^ X ^w, .F}j 



= a,Wj : E, X Fj, (47) 

(-i,i)e/xj 

where ai is the probability of agent Ei (the same for Wj and Fj). 

Expressions ( I46> -( I47> are simple examples of equational laws of WSCCS. A complete development can be found in 
yS.Ul. However, our simple presentation allows to point out the power of WSCCS for society modeling. Hence, let us 
consider the simple example of a colony of ants (agents) that can be only Passive or Active. In this example, described in 
1321 . the active agent is defined by an expression analogous to Equation ( I42> : 

Active = p : ■/ .Passive + {1 — p) : / .Active. (48) 
The passive agent works differently. Following 1^ . we assume that it remains passive forever, thus: 

Passive = 1 : / .Passive. (49) 

The natural question now is: How to combine ants in order to define a colony? This question is answered by the compo- 
sition operation (Expressions ( I46t .( l47» . Henceforth, we write a colony of n ants, i of which are Active agents, as: 

i n—i 

Colonynii) = Active x x Active x Passive x x Passive = Actwe x Y\ Passive (50) 

i agents n-i agents 



Following Expression J47> . we can demonstrate that (see 1321 . page 171, for details): 

k n~k 

k 

fe=0 

We shall obtain the meaning of the coefficients: 



i / . \ / k n~k \ 

Colony,, («) = XI ( I ) P'^'' ('^^Pf ■ ^- [Yl^ctive x J| Passive I 



(51) 



^^>'^ = E I jp'^'^i^-pf- (52) 
Fkstly, according to Equation (I45> . the transitions are given by: 

Colonynii) / [ci.j] Colonyn{j), j = 0, (53) 

> 

In order to interpret Ci,j, we consider now the sequence of random variables A = {At : t £ {0, 1, 2, ...}} where < At < 
n for each t, that represent the outcome of a series of transitions on the agent Colonyn{n), with initial state consisting of all 
ants in the active state. We can think t as the number of ticks of a global clock. From expression i5l\ it is straightforward to 
observe that: 

P{At+i^j\At^i)^c,,j. (54) 



^From this expression, we observe the WSCCS model, given by Equation i5l\ . has underlying discrete time Markov chain. 
A Markov chain is a time ordered sequence of random variables where the t + 1 variable of the sequence is conditional only 
on the tth variable 's value 1 19 1. In fact, this happens for WSCCS models in general (see Appendix A of 1 32 1). Such feature 
is used in 1321 in the context of ant societies. Basically, the transition rules can demonstrate properties that can help the 
analysis of important behaviors (asymptotic ones, for instance). 



6. Lattice Gas Automata and Multiscale Analysis 



The WSCCS is useful for modeling and analysis of the discrete dynamics of agent system. The analysis does not attempt 
to get spatial distribution of observables. Such goal can be achieved by multiscale techniques. In this section we consider 
the FHP model, which is a Lattice Gas Cellular Automata model, used for fluid simulation. Thus, space variables must be 
considered, that means. In this case, a multiscale technique based on Chapman-Enskog |25| expansion is used to establish 
the connection between the microscopic dynamics and the macroscopic observables. 

The Chapman-Enskog method works as follows. Given an operator ^ and the equation: 

e (/) - 0, (55) 

let us suppose that: 

1. The solution / can be expressed as: 

/^/(o)+;(i)+j(2) + ... (56) 

2. When this series is introduced in Expression i55i the result can be expressed as: 

? (/(°) + /(I) + /(^) + ...)= (/(°)) + (/(°\ /(^)) + (/(°\ /(^\ /(^') + ... (57) 



3. The functions are such that: 

C(°)(/(°))=0, (58) 

^(i)(/(o)^y(i))^0, (59) 

C(')(/(°\/«,/('))=0, (60) 



(61) 

which together ensure that Expression i55\ is satisfied. 

Therefore, following items (l)-(3) we say that the sub-series f^-^\ f^^^ + f^^\ /^°-' + /^^■' + f^'^\ are successive 
approximations of /. Arbitrary elements may enter into the solution of Equations as well as in the definition of the 

approximations and of the expansion ( I56> . An interesting example is given by the FHP model. 

The FHP was introduced by Frisch, Hasslacher and Pomeau 1 17 1 in 1986 and is a model of a two-dimensional fluid and it 
is an abstraction, at a microscopic scale, of a fluid. The FHP model describes the motion of particles traveling in a discrete 
space and colliding with each other. The space is discretized in a hexagonal lattice. 

The microdynamics of FHP is given in terms of Boolean variables describing the occupation numbers at each site of the 
lattice and at each time step (i.e. the presence or the absence of a fluid particle). The FHP particles move in discrete time 
steps, with a velocity of constant modulus, pointing along one of the six directions of the lattice. The dynamics is such that no 
more than one particle enters the same site at the same time with the same velocity. This restriction is the exclusion principle; 
it ensures that six Boolean variables at each lattice site are always enough to represent the microdynamics. 

In the absence of collisions, the particles would move in straight lines, along the direction specified by their velocity vector 
The velocity modulus is such that, in a time step, each particle travels one lattice spacing and reaches a nearest-neighbor site. 

In order to conserve the number of particles and the momentum during each interaction, only a few configurations lead 
to a non-trivial collision (i.e. a collision in which the directions of motion have changed). When exactly two particles enter 
the same site with opposite velocities, both of them are deflected by 60 degrees so that the output of the collision is still a 



zero momentum configuration with two particles. When exactly three particles collide with an angle of 120 degrees between 
each other, they bounce back to where they come from (so that the momentum after the collision is zero, as it was before the 
collision). Both two- and three-body collisions are necessary to avoid extra conservation laws. Several variants of the FHP 
model exist in the literature I10II15I . including some with rest particles like models FHP-II and FHP-III. 

For all other configurations no collision occurs and the particles go through as if they were transparent to each other. 

The full microdynamics of the FHP model can be expressed by evolution equations for the occupation numbers defined 
as the number, rii (r, t), of particle entering site r at time t with a velocity pointing along direction Ci, where i = 1, 2, . . . , 6 
labels the six lattice directions. The numbers ni can be or 1. 

We also define the time step as At and the lattice spacing as A^. Thus, the six possible velocities Vi of the particles are 
related to their directions of motion by 

Vi = —Ci. (62) 
Without interactions between particles, the evolution equations for the rii would be given by 

n,{r + ArC„t + At):=:n,{r,t) (63) 

which express that a particle entering site f with velocity along q will continue in a straight line so that, at next time step, 
it will enter site r + ArCi with the same direction of motion. However, due to collisions, a particle can be removed from its 
original direction or another one can be deflected into direction 

For instance, if only m and ni+3 are 1 at site r, a collision occurs and the particle traveUng with velocity Vi will then move 
with either velocity Vi-i or Vi+i, where i ~ 1,2, ... ,6. The quantity 

Di = niUi+s (1 - Ui+i) (1 - ^1+2) (1 - f^^+4) (1 - "1+5) ■ (64) 

indicates, when Di = 1 that such a collision will take place. Therefore ri i — Di is the number of particles left in direction Ci 
due to a two-particle collision along this direction. 

Now, when n,; = 0, a new particle can appear in direction ci, as the result of a collision between ni^i and ni+4 or a 
collision between n^-i e ni^2- It is convenient to introduce a random Boolean variable q (r, <), which decides whether the 
particles are deflected to the right (q — 1) or to the left (q = 0), when a two-body collision takes place. Therefore, the number 
of particle created in direction ci is 

+ (65) 

Particles can also be deflected into (or removed from) direction Ci because of a three-body collision. The quantity which 
express the occurrence of a three-body collision with particles Ui, 71^+2 and ni+4 is 

Ti = nini+2ni+4 (1 - "i+i) (1 - n^+a) (1 - 71^+5) (66) 

As before, the result of a three-body collision is to modify the number of particles in direction ci as 

-T,+ T,+3, (67) 

Thus, according to our collision rules, the microdynamics of a LGCA is written as 

n,{f+ArC„t + At)^n, {r, t) + VL, [n {f, t)) (68) 

where fli is called the collision term. 

For the FHP model, fi^ is defined so as to reproduce the colUsions, that is 

n, = -A + <zA-i + {l-q) A+1 -T, + r,+3. (69) 
Using the full expression for Di and Ti, given by the Equations ( I64> - (I66> . we obtain, 

(70) 

-nini+2»T-i+4 (1 - n^+i) (1 - "-j+3) (1 - rii+5) 
+ nj+inj+3ni+5 (1 - n.,) (1 - ni+2) (1 - n.1+4) 
- (1 - rii+i) (1 - ni+2) (1 - ^^1+4) (1 - ^^^+5) 

+ (l-q) Ui+irii+i (1 - Ui) (1 - ni+2) (1 - "1+3) 
+ (1 - g) (1 - n,+5) 

+ gni+2"j+5 (1 - rii) (1 - "-i+i) (1 - ^i^+a) (1 - rii+i) ■ 



These equations are easy to code in a computer and yield a fast and exact implementation of the model 

Until now, we deal with microscopic quantities. However, the physical quantities of interest are not so much the Boolean 
variables rii but macroscopic quantities or average values, such as, for instance, the average density of particles and the 
average velocity field at each point of the system. Theses quantities are defined from the ensemble average Ni (r, t) — 
{rii (r, t)) of the microscopic occupation variables. Note that, Ni ( f, t) is also the probability of having a particle entering 
the site r, at time t, with velocity 

Vi = -7-Ci- 

In general, a LGCA is characterized by the number z of lattice directions and the spatial dimensionality d. In our case 
d = 2 and z — 6. Following the usual definition of statistical mechanics, the local density of particles is the sum of the 
average number of particles traveling along, each direction Ci 

z 

p(f,i) = ^7V,(f,i). (71) 
Similarly, the particle current, which is the density p times the velocity field u, is expressed by. 

z 

pir,t)u{r,t) = J2v,N,{f,t). (72) 

i=0 

Another quantity which will play an important role in the up coming derivation is the momentum tensor 11 defined as 

z 

na/3 = ^ VtaVipNi [r, t) (73) 
1=0 

where the Greek indices a and /? label the d spatial components of the vectors. The quantity 11 represents the flux of the 
a— component of momentum transported along the /?— axis. This term will contain the pressure contribution and the effects 
of viscosity. 

The starting point to obtain the macroscopic behavior of the CA fluid is to derive an equation for the N'^s. Averaging the 
microdynamics j68> yields 

(f + A.c- , i + At) - (f, t) = {a, {n {r, t))) (74) 

where fli is the collision term of the LGCA, under study. It is important to notice that ^li (n) has some generic properties, 
namely 

z z 

e J2v,a,^0 (75) 

i=l i=l 

expressing the fact that particle number and momentum are conserved during the collision process (the incoming sum of 
mass or momentum equals the outgoing sum). 

The Ni's vary between and 1 and, at a scale L >> A^ e T >> At, one can expect them to be smooth functions of the 
space and time coordinates. Therefore, Equation can be Taylor expanded up to second order and gives 

A, (c, • V) N, (f, t) + AtdtN, (r, t) (76) 

+ i (A,)' (c, • V)' N, (f, t) + ArAt (c, • V) dtN, (f, t) 

+ liAtf (<9t)'iV,(f,<) = (a (n(f,t))). 

where (St)^ is the second derivative in respect to the time parameter t. 

At a macroscopic scale L >> Ar, following the procedure of the so-called multiscale expansion 1281 . we introduce a new 
space variable fi such that 

r*! = edf^ e dr — e.dp^ (77) 

with e << 1. We also introduce the extra time variables ti and t2, as well as new functions depending on ri, ti and t2, 
Nf — {ti,t2,ri) and substitute into Equation (I76> 

N,^N^ dt ^ edt, + e^dt, dr edp, (78) 



together with the corresponding expressions for the second order derivatives. Then obtain new equations for the new functions 
Nf. Thus, following step (1) above we may write [28 1 (see Expression \56V ). 



Nl = N^°^ + eivf ^ + e^TVf ) + • • ■ (79) 

The Chapman-Enskog method is the standard procedure used in statistical mechanics to solve an Equation like ( I76> with 
a perturbation parameter e. Assuming that {Hi (n)) can be factorized into fli (N), we write the contributions of each order 
in e. According to multiscale Expansion j79> . the right-hand side of ( I76> reads 



(80) 



Using Expressions (^)-(29} in the left-hand side of ( 17 6> and comparing the terms of the same order in e in the Equation dSOI l, 
yields 

O (e") : a, (tvW) = (81) 

and 

O (e^) : ^lc.v^c.N^°^ + dt,Nl°^ (82) 



where the subscript 1 in spatial derivatives (e.g. dia) indicates a differential operator expressed in the variable fi and 
{ci ■ Vri) ~ diaVia, from Equation ( I62L 

We also impose the extra conditions that the macroscopic quantities p and pu are entirely given by the zero order of 
Expansion ( I79> 

z z 
1=1 i=l 

and therefore 



^0 and ^ v,N^''^ ^ 0, for ? > 1 (84) 

i=l i=l 

Thus, following the Chapman-Enskog method we can obtain I16II10I . from Equation ( I76> . the following result at order e 

dti p + divi pu — (85) 

and 

dt.pu^ + di^n'^^l^O (86) 
On the other hand, if we considered the terms of order and using the Relations ( I85> and ( I86> to simplify, we have 



dt2PUa + dip 



n(i) i ^* n^'') J- ;) 



(87) 



The last equation contains the dissipative contributions to the Euler Equation ( I86t . The first contribution is Il^^j which is the 

dissipative part of the momentum tensor. The second part, namely {^ti^'^ap + ^17'^q^\) comes from the second order 
terms of the Taylor expansion of the discrete Boltzmann equation. These terms account for the discreteness of the lattice and 
have no counterpart in standard hydrodynamics. As we shall see, they will lead to the so-called lattice viscosity. The order e 
e can be grouped together to give the general equations governing our system. Summing Equations ( I85> and ( I87> with the 
appropriate power of e as factor and we obtain 

dtp + div pu = (88) 



Similarly, Equation ( I86t and ( I87t yields 1161 



dtpUa 



dr 



n, 



Q/3 



At 
2 



= 



(89) 



We now turn to the problem of solving Equation iSU together with conditions i83l in order to find iy-^^-* as functions of p 
and pu. The solutions iV^°'' which make the colUsion term fl vanish are known as the local equilibrium solutions. Physically, 
they correspond to a situation where the rate of each type of collision equilibrates. Since the collision time At is much smaller 
than the macroscopic observation time, it is reasonable to expect, in first approximation that an equilibrium is reached locally. 
Provided that the collision behaves reasonably, it is found 1161 that the generic solution is 



(0) 



1 + cxp ( -A - B ■ V, 



(90) 



This expression has the form of a Fermi-Dirac distribution. This is a consequence of the exclusion principle we have imposed 
in the cellular automata rule (no more than one particle per site and direction). This form is explicitly obtained for the FHP 
model by assuming that the rate of direct and inverse collisions are equal. The quantities A e i? are functions of the density 
p and the velocity field u and are to be determined according to Equations ( I83> . In order to carry out this calculation, N^^'' is 
Taylor expanded up to second order in the velocity field u. One obtains 1101 



aAo) bp^ , pG{p) 

N> ' = ap+ —Vi ■ U H -^QiapUaU/S 



(91) 



where a, f3, 7 are summed over the spacial coordinates, e.g. a, (3, j & {1, ... , d}, v 



-, a 



i,6 



and 



-5a 



(92) 



The function G is obtained from the fact that Nf'^ is the Taylor expansion of a Fermi-Dirac distribution. For FHP, it is 
found iTOIITel 

2 (3 - p) 



G{p) = 



3 (6 - p) 



We may now compute the local equilibrium part of the momentum tensor, li'^^l and then obtain the pressure term 



p = aC2v'^p 

where C2 = f . 

We can see 1161 that the lattice viscosity is given by 



C2 „ 



pGip)u' 



(93) 



^lattice 



-C46 



Atv^ 



dA 



d{d + 2) z 2 
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V 



-A 
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V 



2{d + 2) 

The usual contribution to viscosity is due to the collision between the fluid particles and is given by fW\ 



where —A is given by —A = 2s (1 — s)'^ where s = | 
Therefore, the Navier-Stokes equation reads 

dtu + 2CiGip) (u- V)m 



1, 

P 



-Vp + v^^u 



(94) 



where 













(1- 




d + 2 




-I) 



(95) 



is the kinematic viscosity of our discrete fluid. 

Therefore, we demonstrated that the Navier-Stokes model can be reproduced by FHP technique. However, there is no 
need to solve Partial Differential Equations (PDEs) to obtain a high level of description. Such advantage can be explored in 
technological and scientific applications. For instance, in 1181 we propose to combine the advantage of the low computational 
cost of LGCA and its ability to mimic the realistic fluid dynamics to develop a new animating framework for computer 
graphics applications. 

7. Tool for Agent-Based Simulation 

Agent-based models can be analyzed by computer simulations. The NetLogo software is one possibility in this area 
1271 . It is a programmable modeling environment for simulating complex systems developing over time. Modelers can give 
instructions to hundreds or thousands of independent " agents" all operating concurrently in order to explore the connection 
between the behavior of individuals and the macroscopic patterns that emerge from the interaction of many individuals. 
Users can create their own models using NetLogo facilities and documentations. It also comes with a Library of pre- written 
simulations that can be used and modified. 

As an example of the NetLogo capabilities we describe our implementation of a Lattice Gas model called HPP II II . It is 
similar to the FHP model described on Section|6]but, in this case, the lattice is a rectangular one. Figure|5]shows the NetLogo 
main interface and a snapshot of our HPP implementation. 



NetLogo: example(HPP)_01-mod {Z:\netlogo-2.1} 
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Interface Plnlormatior|| Procedures]! EttQ's| 
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Plot S3 Output 



max valuesTrorrr patches [count turtles-here] 
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Figure 5. The NetLogo main interface and our l-IPP implementation 

The rules used for collision are explained in Figure |6l In the other situations the particles are considered transparent to 
each other when they cross the same site. There is also an exclusion principle: it is not allowed more than 1 particle entering 
a given site with a given direction of motion. The aim of these rules is to reproduce some aspect of the real interactions 
between particles, namely that momentum and particle number are conserved during a collision. With such simple dynamics, 
we can model and simulate a gas of colliding particles and to obtain complex behaviors [ 11 1. 

The HPP model is a kind of cellular automaton which has a lattice of sites that may have 0, 1, 2, 3 or 4 crossing particles 
at a time t. The rules define the system (particles) evolution and, consequently, the update of each site value. 

The evolution of the sites is often split in two steps: collision and motion (or propagation). The collision phase solves 
interactions (collisions) through the rules pictured on Figure|6| During the propagation phase, the particles actually move to 
the nearest neighbor site they are traveling to. 



(a) 



Ik 

(c) 

Figure 6. HPP rules for collision. 

The implementation the HPP in the NetLogo software we must define the agents, which are represented by arrows in the 
Figure|6] and the lattice. In the NetLogo system, the " bricks" to compose an application are: 

1 . Application Control: Button, Slider, Switch, Chooser. 

2. Plot. 

3. Monitor, Output and Text. 

4. Turtles: agents plus their graphical representation. 

The Figure Is] shows the instances of some of these tools in our HPP implementation. The basic controls for the model are 
the following Buttons: (a) SETUP - Sets up screen with a given percentage of particles; (b) Execute - Run the model; (c) 
MOVE TURTLES- For move the particles with the mouse. There is one Slider to set the number of particles. 

Behind the graphical interface for visualization and control the application, there is a code that implements agents behav- 
iors. For example, let us consider the following code line: 

if any! other — arrows — here with (96) 

[heading = heading — of myself and who < who — of myself] [jump — 1] 

if: Reporter must report a boolean (true or false) value. 

any?: Reports true if the given agentset is non-empty, false otherwise. 

other-BREED-here: Reports an agentset consisting of all turtles on the calling turtle's patch (not including the caller 
itself). If a breed (a built-in turtle variable) is specified, like arrows in the above example, only turtles with the given breed 
are included. 

heading: It is command in the NetLogo syntax. Each turtle picks a random integer between and 359. Then the turtle 
sets its heading to the number it picked. Heading is measured in degrees, clockwise around the circle, starting with degrees 
at twelve o'clock (north). 

myself: It means " the turtle or patch who asked me to do what I'm doing right now". 

wlio?: This is a built-in turtle variable. It holds the turtle's id number (an integer greater than or equal to zero). You 
cannot set this variable; a turtle's id number never changes. When NetLogo starts, or after you use the clear-all or clear- 
turtles commands, new turtles are created with ids in order, starting at 0. If a turtle dies, though, a new turtle may eventually 
be assigned the same id number that was used by the dead turtle. 




jump: This is another command. Turtles move forward by number units all at once, without the amount of time passing 
depending on the distance. 

NetLogo system has a lot of examples and a good documentation to help new users to write its own applications. 

8. Conclusions 

The simulation of dynamical systems through a large set of interacting agents is an interesting research field with appUca- 
tions in areas like, physics, economy and sociology. This is a bottom-ut approach which tries to derive global properties of a 
complex system through local interaction rules and agent behavior. Agent-Based Modeling has the advantage of simphcity 
and low computational cost if compared with the traditional differential equation approaches. 

In this paper we survey a method based on the WSCCS to express agent behavior which allow to analytically predict 
the results obtained by the simulation. Also, multiscale techniques, based on Chapman-Enskog expansion was reviewed 
to establish the connection between the microscopic dynamics (agent behavior) and the macroscopic observables. Besides, 
Principal Component Analysis (PCA) was analyzed for knowledge discovery in a Cellular Automata database. Finally, we 
show the capabilities of the NetLogo, a free software for agent simulation of complex system and describe our experience 
with this package. Our research will continue in this field, specially exploring the appUcation of agent-based models for 
computer graphics applications. 
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